QCD Evolution by Finite Element Methods 
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Abstract 

A simple, new method for solving for the Q^ evolution of parton 
distributions in perturbative QCD using cubic splines is described and 
applied to the evolution of nonsinglet quark distributions. 

1 Introduction 

The efficacy of using perturbative QCD to describe processes at large mo- 
mentum transfer in terms of process independent parton distributions has 
spawned an industry devoted to the determination of these distributions |I|]. 
Central to this endeavour is the ability to solve the Altarelli-Parisi equations 
which govern the Q^ evolution of parton distributions, so that data taken at 
different energies can be used to constrain the form of the parton distribu- 
tions. The evolution equations have been solved in a variety of ways, ranging 
from brute force integration of the integro-differential equations^, to orthog- 
onal polynomial expansions 0], to solutions in terms of Bernstine polynomials 
for next to leading order evolution^. In this report, a straightforward new 
method for solving the Altarelli-Parisi equations in terms of a small number 
of cubic splines is described and illustrated for the case of nonsinglet quark 
distributions. Since the method relies only on the smoothness of the parton 
distributions and the moment structure of the QCD evolution equations, it 
is applicable to singlet and nonsinglet distributions in both leading and next 
to leading order with equal ease. 



2 The Method 

In this section, the apphcation of finite element methods to the problem 
of QCD evolution is demonstrated for the simple case of a nonsinglet quark 
distribution. Given a nonsinglet quark distribution qNs{x, Ql) at some initial 
renormalization scale Qq, we are looking for a new distribution evolved to a 
different Q^ according to QCD perturbation theory. To leading order, this 
requires 

-M„(g^) = -M„(g^)x(|^)'^"^ (1) 

where AiniQ"^) = Jo dxx"-"^ (1ns{x,Q'^) is the nth moment of the quark 
distribution, and d^^ is the anomalous dimension of the nth nonsinglet, twist 
two operator. 

In principle, knowledge of all the moments is required to uniquely de- 
termine qNs{x,Q^), but the distributions used in QCD phenomenology are 
peaked at small values of x, and all but a few of the moments are small. Ex- 
ploiting this fact, we proceed by numerically calculating the first N moments 
of the desired distribution in terms of a set of N suitably chosen integration 
points and weights. Thus, 

N 

MniQ') ^J2^r'wms{x^,Q'), (2) 

where {xi} are the integration points and {wi} are the weights. Assuming 
that the integrals are well approximated, this yields a matrix equation relat- 
ing the first A^ moments of the distribution to the values of the distribution 



function at the points {xj}. In principle, the distribution can be solved for 
by taking the inverse of the matrix and the limit as the number of integra- 
tion points gets large. As a practical method, this procedure fails since the 
matrix, 

is difficult to invert numerically when N is large. For relatively small values 
of A^, however, Ani can be inverted easily and quickly by any one of a variety 
of standard techniques. The result is an approximation of the values of the 
parton distribution at A^ points. All that is left is to fill in the rest of the 
curve between these points. Since the distribution is expected to be smooth 
and non-oscillatory, the remainder of the curve is found by interpolating 
between the points Xi using cubic splines[^]. Since the functions of interest 
are singular at x = 0, some care must be taken when choosing the integration 
scheme and spline basis so that the function may be adequately fit at small x. 
While other techniques for dealing with this problem exist in the literature0], 
here we simply change the variable in the numerical integration and spline 
fit to y = x", with a < 1, so that small values of x have increased weight in 
the integrations and the curve is better approximated by a polynomial in y 
near x = 0. 



3 Results and Discussion 

The essence of the method outhned in the previous section is the reconstruc- 
tion of an essentially arbitrary curve from its moments. The simplest test 
available is to decompose a curve into moments and then reconstruct it with- 
out any change in Q^ . In Figure 1, the results are displayed for a curve typical 
of those encountered in QCD phenomenology, xqNs{x) = y/x{l~x)^, using a 
Gauss-Legendre integration scheme with a = 2/3, N=3,5,7 or 10, and splines 
that are functions oi y = ^/x. For N=3 the procedure has clearly failed to 
reproduce the original, while for N=5 the reconstructed curve deviates only 
slightly. For N=7 and 10, however, the reconstructed curves are virtually 
indistinguishable from the original. In Figure 2, a set of 10 splines is used 
to perform the leading order evolution of the valence d quark distribution 
parametrized by Morfin and Tung(MT)[^] from an initial scale Qg = 4 GeV^ 
to Q^ = 10 and 50 GeV^. Using a = 1/2 for both the integration points and 
splines, the MT curves are well reproduced. (The small deviation from the 
MT curve is a reflection of the approximate nature of the parametrization of 
the Q^ dependence of xdv{x) in reference 4.) 

The finite element method for solving the Altarelli-Parisi equations has 
a number of conceptual and practical advantages. To begin, the method is 
conceptually straightforward, relying only on the moment structure of the 
evolution equations and the smoothness of the quark distributions. Since 
these properties hold for all parton distributions the method is applicable, 



essentially without change, to both singlet and nonsinglet distributions in 
both leading and next to leading order in as- In addition, one need not com- 
mit to a particular form for the parton distributions when fitting data, but 
can instead treat the moments of the distribution, which are related to the 
physically relevant twist two matrix elements, as fit parameters. Practically, 
the spline methods used here are more stable than schemes involving orthogo- 
nal polynomials, which are prone to global oscillations when singularities are 
encountered. Since the elements of the scheme we have described, numerical 
integration, matrix inversion, and spline interpolation, are extensively used in 
other applications, the method is easily implemented. If greater accuracy is 
desired, the method can be altered to include any non-integral, finite moment 
of the parton distributions by simply generating the appropriate anomalous 
dimensions and coefficients. Since the matrix Ani need only be inverted once 
to fit data over a large range in Q^ , the current method is faster than brute 
force integration of the Altarelli-Parisi equations. 
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Figure Captions 

Figure 1 Reconstruction of the curve ^/x{l — xY using 3,5,7 and 10 
splines. 

Figure 2 Leading order evolution of the valence d quark distribution of 
Morfin and Tung(MT)||]. 
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